Exact treatment of planar two-electron quantum dots: effects of anharmonicity on the 
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Static properties of an anharmonic potential model for planar two-electron quantum dots are 
investigated using a method which allows for the exact representation of the matrix elements, in- 
cluding the full Coulombic electron - electron interaction. The anharmonic confining potential in 
combination with the interparticle Coulomb interaction affects the spectral properties of the sys- 
tem considerably as it implies total loss of separability of the system. Properties of the classical 
phase space, spectral measures of the chaoticity, as well as localization properties of the eigenstates 
corroborate this. 
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I. INTRODUCTION 

The understanding of the complexity of a physical 
system has been addressed for the last decades and its 
physical relevance is still a matter of investigation. For 
open systems the correlation of quantum chaos and trans- 
port is widely discussed, e.g. for quantum dots 1 , giving 
also evidence for the onset of decoherence and the quan- 
tum to classical transition 2 . The connection between 
thermalization and classically chaotic dynamics is well 
known since the investigations of ergodicity in the Sinai 
billiard 3 , while the link between quantum thermalization 
and quantum chaos is not clear yet 4 . Experimental tech- 
niques for the coupling of microcavities to waveguides 
exploit the underlying chaoticity of the system . 

One of the simplest realizations of a complex system in 
atomic physics is the three-body Coulomb problem. The 
electron-electron interaction term in the Hamiltonian of 
the unperturbed helium atom renders the two-electron 
dynamics in general chaotic with only rather small re- 
gions of the classical phase space occupied by regular 
motion. On the quantum level, the loss of integrabil- 
ity leads to an abundance of intriguing and surprising 
effects 6-8 . Replacing the attractive Coulomb potential 
by an alternative kind of potential generates different, 
though in many cases equally challenging physical mod- 
els. 

The Hooke's atom - that is two electrons with har- 
monic confinement - has been thoroughly investigated, 
not only because it serves as the most common model 
for a generic semiconductor quantum dot. Its relative 
motion can be separated from the center of mass (COM) 
motion, a fact generally known as Kohn's theorem 9 . Nev- 
ertheless, the classical dynamics of the three-dimensional 
system is in general mixed regular-chaotic 10 . On the 
quantum level, there is a reduced family of eigenstates, 
the energies of which admit closed analytical solutions 
for special values of the confining harmonic potential, 



the Coulomb interaction between the particles and pos- 
sibly an additional magnetic field perpendicular to the 
plane 11 . 

Due to its simplicity, the harmonic two-electron quan- 
tum dot has been used as a paradigmatic model for 
the study, for example, of the entanglement of two 
electrons 12,13 or of the origin of Hund's rule 14 . In the 
former case, the entanglement of low-energetic states in- 
creases with the interaction between the electrons and 
with the energy of the state. This is consistent with in- 
vestigations in low-excited states of helium 15 . However, 
the onset of chaotic dynamics might induce qualitatively 
new features 16 . In the latter case, the spin symmetry de- 
termines the symmetry of the spatial wave function un- 
der particle exchange, leading to the effect known as the 
Fermi hole, describing the minimum of the triplet wave 
function around the origin. Still, it is not a priori clear, 
which cigenenergy is energetically favorable for identical 
spatial configuration but different spin symmetry 

Extending the harmonic model, the next to leading or- 
der quartic term in the potential has so far received very 
little attention in the literature, e.g. Ref. 17 and 18. 
Yet, the enhanced complexity of the dynamics induced 
by this anharmonicity - as a consequence of the loss of 
separability between the COM and the relative motion - 
might help for a better understanding of, e.g., the entan- 
glement of two electrons in atomic systems 19 . Further- 
more the validity of Hund's rules in an anharmonic case 
has not been addressed, so far. A similar effect has been 
seen earlier for the quantum dot model, where the spin 
symmetry of the ground state depends on the applied 
magnetic field 20 . 

Measurements of the electric current through a single 
quantum dot depending on the applied gate voltage show 
a specific shell structure of the ground state energies for 
few electrons confined in the dot 21 ' 22 . This shell struc- 
ture supports a planar approach with a harmonic con- 
fining potential, which has extensively been studied in 
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the past, most frequently for the two electron case 23 . A 
detailed review on the properties and modeling of semi- 
conductor quantum dots can be found in Ref. 24. The 
question of the relevance of the planar model has been 
discussed in detail in Ref. 10, including the effects of 
a static magnetic field perpendicular to the plane. The 
planar confinement generates an energy shift, which can 
be related to the ground state energy of a strongly con- 
fined harmonic potential in the z-direction, E z = ftio z /2. 
In addition, the Coulomb interaction is overestimated by 
the planar restriction. These effects can be properly in- 
corporated in a planar approach by an appropriate rescal- 
ing justifying this model under strong confinement to the 
plane, uj z 3> co xy . 

The aim of this work is to shed some light on the under- 
standing of the complex dynamics induced by a quartic 
potential in a two-electron quantum dot. For this pur- 
pose we have developed a quantum approach to an an- 
harmonic two-electron quantum dot confined to a plane 
in analogy to a planar model for two-electron atoms 25 ' 26 . 
In Sec. II we describe our extended model of a two- 
electron quantum dot and clarify the underlying sym- 
metry structure, already inherent to the non-interacting 
harmonic model. We show the convergence of our results 
and compare to available data in the literature. After a 
brief description of the mixed regular-chaotic structure 
of the classical phase space of the system in Sec. Ill A, 
we proceed to explore the role of the anharmonicity in 
the complexity of the quantum mechanical model in Sec. 
IIIB. Finally we conclude in Sec. IV. 



II. THEORETICAL APPROACH 

A. Physical model 

The problem we wish to describe with a minimum of 
approximations is a system of two interacting fermionic 
negative charges confined in an anharmonic potential. 
Each charge has an effective mass m* and the Coulomb 
repulsion between them is affected by the dielectric prop- 
erties of the confining system manifested through the di- 
electric constant e. Introducing modified atomic units 
(see Appendix A) the quantities m* and 47reoe can be 
set to unity. In the course of this work we are going to 
use these modified units ([a.u.]). The general form of the 
Hamiltonian describing our model reads as 



H = J2 



-iv 2 

2 J 



+ V con {(rj) 



+ V int (r 12 ). (1) 



the dynamics in the quantum dot is strongly confined to 
the xy '-plane. This is the case we consider in this contri- 
bution. From now on, we restrict the dynamics to two 
dimensions of configuration space, with the Cartesian po- 
sitions (xi,yi) and (22,2/2) of the charges. We assume 
furthermore that the confinement potential is isotropic, 
V con {(r) = \lo 2 t 2 + Kr 4 (cj = uj x = u v ). 

The eigenvalue problem for 7 7^ effectively depends 
on two parameters, namely k and ui. Following experi- 
mental data we set uj = 1 (see appendix A). For small 
values of the quartic potential strength k the experimen- 
tally confirmed symmetry structure of the planar har- 
monic model is effectively not broken 21 . Therefore we 
vary the anharmonicity of the potential with the help of 
the parameter k, which ranges from 0.0 to 0.1 for the 
interacting (7 = 1) and non-interacting cases (7 = 0). 



B. Underlying symmetry structure 

To gain insight into the level structure of the full sys- 
tem (1), it is instructive to consider the energy levels 
and their degeneracies in the fermionic, non-interacting 
planar harmonic oscillator model. There are several pos- 
sibilities to distinguish between symmetry classes, but 
our choice is made in order to resolve all degeneracies as 
soon as any kind of perturbation is introduced (re ^ or 
7^0). 

We start from the radial representation in COM and 
relative coordinates, which gives rise to the eigenbasis 
\n c , m c , n r , m r ) with the principal quantum numbers n c 
and n r and the angular momentum quantum numbers 
perpendicular to the plane m c and m r , for the COM and 
relative motion, respectively. The energies are given by 
E c i r = u)(2n c i r + \m c / r \ + 1). These quantum numbers 
lose their meaning when the anharmonic potential is in- 
troduced due to the coupling of the COM and relative 
motion. On the other hand, the squared total angular 
momentum operator perpendicular to the plane L 2 , the 
particle exchange operator II12 and the operator IL xy , 
which interchanges the spatial coordinates Xi and j/i, 
commute with the Hamiltonian (1). The latter is a two- 
dimensional parity operator, as it changes the orientation 
of the coordinate system. The action of the exchange and 
interchange operators on some function ip(xi,y±,X2,y2) 
in coordinate space is given by 

II12 2/i> 2=2, 2/2) = ±ip{x 2 ,y 2 ,x 1 ,y 1 ), (2) 
Kx y ip(zi,yi,X2,V2) = ±tp{yi,x 1 ,y 2 ,x 2 ), (3) 



Here Vj acts only on rj and Vj n t (1*12) = 1/ r i2 is the inter- 
action potential between the two charges separated by a 
distance r\ 2 = |ti — r 2 \- For practical purposes we intro- 
duce a control parameter 7 such that Vi n t(?"i2) = 7/V12. 
The anharmonicity of the quantum dot is modeled by 
a quartic contribution Kr 4 to the harmonic confinement 
iw 2 a; 2 + k^ly 2 + h^l 2 " 2 - If we assume oo z ^> tj x ,ui y , 



respectively. From this and from the coordinate repre- 
sentation of L 2 , 



{xiP Vl - yip Xl + x 2 p V2 - yzP X2 Y 



(4) 



it becomes apparent, that these three and the full Hamil- 
tonian are mutually commuting operators. A common 
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set of eigenstates of L\ , II12 and Tl^y is defined by 

\n c ,m c ,n r ,m r ) e *> = 

(\n c ,m c ,n r ,m r ) + e p \n c -m c n r -m r ))/ 'v2, 

with e p e {+1 (even), —1 (odd)}. In order to guarantee 
uniqueness of the representation of this eigenbasis it is 
necessary to impose m c > and if m c = then m r > 0. 
For the case with m c = m r = 0, we set: 



n c ,m c ,n r ,m r 



n c ,m c ,n r ,m r 



The action of the symmetry operators on this basis is 
given by 

L 2 z \n c ,m c ,n r ,m r Y 11 ' = \m c + m r \ 2 \n c ,m c ,n r ,m r Y P , 



U 12 \n c ,m c ,n r ,m r ) e " 
U xy \n Cl m Cl n ri m r ) £p 



\n c ,m c ,n r ,m r ) ep , 
\n c ,m Cl n ri m r yv , 



where the second identity stems from the fact, that the 
particle interchange does only affect the relative coordi- 
nates since it introduces a rotation by tt about the rela- 
tive z-axes. We identify the quantum numbers: 

\m c + m r \ = m £ No, 
2(m r mod2) = s 6 {0 (singlet), 2 (triplet)}, 

where the choice of s for the particle exchange opera- 
tor is motivated by the numerical basis representation 
(see Sec. II D). The principal quantum number n de- 
scribing the unperturbed energy levels E n = uj(n + 2) 
is n = 2(n c + n r ) + \rn c \ + \m r \. The structure of the 
unperturbed spectrum for the lowest lying states of total 
angular momentum and 1 is shown in Table I. 

Turning on the intcrparticle interaction does not intro- 
duce any change for the COM motion. For the relative 
motion there is not any more a closed expression for the 
energy. Nevertheless two quantum numbers h r and rh r 
can be identified using, for instance, a WKB approach 2 '. 
Therefore, the classification is still exact as long as no 
anharmonic interaction is present. 

For the full potential case the principal quantum num- 
ber loses its meaning when avoided crossings appear in 
the spectrum and symmetries of energy eigenstates inter- 
change (see Sec. IIIB2). 



C. Operator representation 

All relevant physical information is contained in the 
spectrum of the Hamiltonian (1). One of the major 
complications for the diagonalization of this Hamilto- 
nian, inherent to all numerical approaches considering 
few-body problems, is the treatment of Coulomb singu- 
larities. These can be rigorously regularized using a rep- 
resentation in the parabolic coordinates (/i+, v+, /i-, v ) 
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TABLE I. The degeneracies of the planar two-dimensional 
harmonic oscillator with two Fermions with respect to the 
symmetry operators L 2 Z , II12 and ITa,,,. The total degeneracy 
H^En = |(n+l)(n + 2)(n + 3) takes into account all symmetry 
classes, whereas K n only counts the degeneracy within the 
specified symmetry class. 



defined by 



/!+ = y/R+ + X+, X± 

v + = sgn(y+) v / ^+ - x+, y± 

R± = + 



X\ ± X2, 

yi ±2/2, 



/i_ = R- + xl , 

V- = sgn(y_)v/i?_ - 



(5) 



R. 



where ^/g is the Jacobian of the transformation. Notice 
that ri2 = R- = A* 2 - + v_ is a polynomial expression 
of the new coordinates. Furthermore, the kinetic energy 
K = — — after multiplication by the Jacobian, 



1 " 

4 . 

0*5 



)(^ + +0 



(6) 



is a polynomial expression of the parabolic coordinates 
and their derivatives. The same holds for all terms of the 
generalized eigenvalue problem (GEVP) 



EBm 



(7) 



with A = ^/gH and B = y 7 ^, obtained after multiplica- 
tion of the stationary Schrodinger equation by the Jaco- 
bian yjg. This offers the opportunity of a representation 
in circular harmonic oscillator creation and annihilation 



4 



operators defined by 
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After normal ordering we arrive at a representation of A 
and B consisting of 2088 and 25 ordered monomials, re- 
spectively, and the maximum degree is 12. For example, 
we give a few terms of the regularized quartic potential: 

iW = \ + ^(4) 3 «i + -^^{alfazal + ... . (8) 

A much simpler expression is obtained for the angular 
momentum L z = x ± p yi - y 1 p Xl + x 2 p V2 - y2Px 2 , 

L * = ^(alai-4 fl 2+ a 3 a 3-al a 4) = -(ni-n 2 +n 3 -n 4 ) , 

where fik = a\a,k are the corresponding number opera- 
tors. 



D. Basis representation 

Since the circular operators satisfy the usual commu- 
tation relations, 

[a l ,a j }=0, [aj,at]=0, [en , at] = Sij , (9) 

for i,j = 1,2, 3, 4, we can associate a harmonic oscillator 
with each pair of circular operators a\ and a%, which in- 
duces a natural basis set composed of tensor products of 
harmonic oscillator Fock states: 

|7ii7i 2 n 3 n4) = |m) (8 \n 2 ) ® \n 3 ) g) |n 4 ). (10) 

Each monomial clement of a polynomial operator O in 
ladder operator representation couples each basis element 
with exactly one element of the basis. Two elements 
\n1n2n3n4) and \ri! 1 n' 2 n' z ri'^ of the basis set (10) are cou- 
pled or satisfy the selection rule {Ani, An 2 , An 3 , A714}, 
with Arii = n, — n^, if (nin2n3n4\0\n' 1 n 2 n 3 n' 4 ) ^ 0. 
For example, the selection rule defined by the mono- 
mial aifa^) 3 0,30,4, appearing in the representation of 
the quartic potential (8), is {Ani, An 2 , A713, A114} = 



{ — 1,3,-1,-5}. The operator A in (7) defines 171 se- 
lection rules, while the Jacobian operator B has 9 selec- 
tion rules. The only selection rule of the angular mo- 
mentum L z is {0,0,0,0} which trivially implies that the 
basis elements \n1n2n3n4) are eigenvectors of L z (with 
eigenvalue \{n\ — n 2 + n 3 — 714)). For a given selection 
rule An = {Ani, An 2 , Ati 3 , Ari4}, the matrix elements 
(n + An\A\n) and (n + An\B\n), with \n) — |nin 2 n 3 7i4) 
and |n + An) = \ni+Ani 7i 2 + At7 2 7i 3 + At7 3 714 + A774 ), 
involve square roots of integer numbers and depend only 
on rii, n-2, n.3 and 71.4. For example, the matrix el- 
ement of the operator A for the selection rule An = 
{ — 1, 3, —1, —5} reads 

(n + An\A\n) = yjm{n 2 + l)(n 2 + 2)(n 2 + 3) 

x \J nz(rii — 3) (774 — 2) (714 — 1)77,4 

( 5 5 1 5 \ 
x i (m - 1) H n 2 ■ 

V 128 256 v ; 512 J 

The parabolic transformation introduces a four-times- 
layer of the coordinate space resulting in unphysical sym- 
metries. This must be compensated by a restriction of 
the allowed basis vectors. Indeed, only even values of 
ni — n 2 and n 3 — 714 have a physical meaning. 

Since a particle exchange can be identified with a ro- 
tation in the parabolic coordinate subspace (/i-, V-) the 
basis (10) is an eigenbasis of Ili 2 . This follows from the 
identity 

U 12 {ii + ,v + ,ii-,vJ) =(/a + ,i/ + ,±z/_,=f/z_) 

with L_ = — i (yU_<9„_ — ^_<9 M _) = a 3 a 3 — a\a^ 7 which 
equivalently can be written as 

IIi 2 I nin 2 n 3 n 4 ) =e ±1 ^ L ~ | nin 2 n 3 n 4 ) 

=e 2 I nin 2 n 3 n 4 ) . 

The two symmetry classes defined by Ili 2 are thus iden- 
tified with the quantum number s = (singlet states) or 
s = 2 (triplet states) such that n 3 — 774 = s (mod 4). 

The action of the coordinate exchange operator H xy 
on an element of the basis (10) is given by 

n 2!y |nin 2 77 3 774) = |n 2 77in 4 n 3 ) . (11) 

This can be easily seen from the coordinate representa- 
tion of the operators hk, k = 1, 2, 3, 4: 

fik = \{ R ±- #±V± - 1) + K-l) fc {x±d y± - y ± d x± ) , 

with V± = d x± + dy ± . Plus signs correspond to k £ 
{1,2}, while minus signs correspond to k e {3,4}. 

Summarizing, a common eigenbasis of the operators 
L 2 Z , II12 and Tl xy is defined by 

|nin 2 n 3 n 4 ) ep = (|nin 2 77 3 7i 4 ) + ep|77 2 rii77 4 77 3 )) , (12) 
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and the associated quantum numbers are 
1 . 

m = - \m - n 2 + n 3 - ra 4 | , 

s = (713 — n 4 ) (mod 4), 
e p = ±1 , 

respectively. 

To ensure a unique representation of the basis vectors 
we impose restrictions on the basis set resulting in three 
cases: 

1) n3 > n 4 ; ni,7i2 arbitrary, 

2) ri3 = ri4 and ni > 7x2, 

3) 713 = n 4 and ni = ri2, 

but then \n\n2n^n^ tv= = \n1n2n3n4) . 



E. Observables and expectation values 

The expectation value (O) of an observable O for a gen- 
eral state is evaluated with the help of the expression 
(^|BO|^). The Jacobian matrix B must be included due 
to the orthogonality relation for the eigenstates of 



the GEVP (7), 



Knowing the expansion 



coefficients of in the symmetrized basis \n} £p (Eq. 
(12)) the expectation value (O) can be readily obtained 
from the matrix representation of BO in this basis. This 
is a simple task in the case that BO is a polynomial 
function of the parabolic coordinates and, therefore, has 
a finite representation in creation and annihilation oper- 
ators. Operators satisfying such property are the relative 
distance r re \ and the COM position i?coM- For example, 



f 



f 



BRcom — 2 B ^+ — iq 



is clearly a polynomial expression of the parabolic coor- 
dinates. Its ordered representation in ladder operators 
contains 70 terms and defines 15 selection rules. The 
cosine of the angle ipn between the electron radii can 
be estimated from a combination of radial expectation 
values: 



COSiy9i2 



(ri • r 2 ) 



(R 2 -) 



The structure of our code allows for separate calculation 
of the different parts of the Hamiltonian, the kinetic term 
T, and the three potential terms involved. It is easy to 
show, that our system satisfies a generalized virial theo- 
rem, 

) + 4(V r quart ic) — (Vcoulomb) , 

where (•) denotes the expectation value in an eigenstate. 
Equivalently this may be changed to an expression for 



the energy depending only on the potential terms of the 
system: 



E = 2(V h 



armonic/ 



3(F C 



1 



quart ic/ 



0'( 



Coulomb / 



(13) 



Our numerical results agree with the virial theorem (13) 
up to the full accuracy of the eigenenergies. 



F. Convergence of the method 

The symmetrized basis (12) allows an analytic, exact 
representation of the GEVP (7). For practical purposes 
this basis has to be truncated. We use the truncation 
criterion 



"1+^2+ m + ni < rib a 



(14) 



for a given integer number nbase- Typical values of nbase 
in our calculations go up to 140, which translate into ba- 
sis dimensions up to 16530. A Krylov subspace method 
for symmetric matrices - which is generally known as 
the Lanczos algorithm 28 - exploits the band structure of 
these matrices and efficiently calculates the largest eigen- 
values of the problem. We take advantage of this prop- 
erty by shifting the energy and solving the inverse prob- 
lem, which leads to a well converged spectrum around 
the shifted energy 29 ' 30 . 

In order to verify the convergence of the eigenvalues ob- 
tained we introduce the scaling transformation mediated 
by the unitary operator P a = cxp[— ^ (r ■ p + p ■ r ) log a] , 
with a a real scaling factor. Position and momentum 
are transformed according to r — » ar and p — > p/a. 
The physical properties are not altered by this unitary 
transformation. In particular the exact eigenvalues are 
invariant under this transformation, 



dE a 
da 



{E a \^\E a 
da 



0. 



(15) 



The truncation of the basis leads however to an a- 
dependency of eigenenergies. The expectation value 
(15) does not vanish anymore, but it is small for well- 
convcrged eigenvalues, a can thus be treated as a vari- 
ational parameter. For convergence of an eigenvalue E a 
we demand that (E a \dH/da\E a ) < 0.1. With this con- 
dition we typically obtain at least three figures of the 
eigenvalues converged. 

This is illustrated in the following with the interacting 
harmonic case n — and with the full potential case for 
ui = 1, k = 0.1 and 7 = 1. 



1. Harmonic case 

There exist analytical solutions for the case without 
quartic potential (to — 1, 7 = 1, k = 0) 11 , which in 
general only appear for very special combinations of the 



6 



E 


(dH/da) 


(-Room) (numerical) 




{Room) (an 


alytical) 


{rid) (numerical) 


2 


999 


999 


999 


999 


999 





000 


000 


000 


000 


049 





626 


657 


068 


657 


750 





626 


657 


068 


657 


750 


1 


636 


801 


341 


900 


274 


5 


000 


000 


000 


000 


005 


-0 


000 


000 


000 


000 


008 


1 


096 


649 


870 


151 


057 


1 


096 


649 


870 


151 


063 


1 


636 


801 


341 


900 


264 


7 


000 


000 


000 


000 


039 


-0 


000 


000 


000 


000 


012 


1 


419 


769 


921 


177 


681 


1 


419 


769 


921 


177 


715 


1 


636 


801 


341 


900 


233 


9 


000 


000 


000 


000 


684 


-0 


000 


000 


000 


000 


206 


1 


681 


692 


992 


842 


792 


1 


681 


692 


992 


843 


259 


1 


636 


801 


341 


899 


816 


11 


000 


000 


000 


001 


018 





000 


000 


000 


000 


223 


1 


907 


854 


079 


694 


065 


1 


907 


854 


079 


694 


655 


1 


636 


801 


341 


899 


770 


13 


000 


000 


000 


000 


819 





000 


000 


000 


000 


113 


2 


109 


832 


797 


669 


925 


2 


109 


832 


797 


670 


389 


1 


636 


801 


341 


899 


910 


15 


000 


000 


000 


000 


076 


-0 


000 


000 


000 


000 


100 


2 


294 


050 


646 


216 


108 


2 


294 


050 


646 


216 


114 


1 


636 


801 


341 


900 


266 


17 


000 


000 


000 


000 


028 





000 


000 


000 


000 


184 


2 


464 


507 


391 


773 


477 


2 


464 


507 


391 


773 


425 


1 


636 


801 


341 


900 


298 


19 


000 


000 


000 


000 


298 





000 


000 


000 


050 


815 


2 


623 


894 


283 


580 


835 


2 


623 


894 


283 


588 


085 


1 


636 


801 


341 


895 


081 


21 


000 


000 


000 


132 


665 


-0 


000 


000 


022 


123 


264 


2 


774 


124 


916 


536 


298 


2 


774 


124 


922 


956 


349 


1 


636 


801 


338 


158 


184 


23 


000 


000 


005 


786 


998 





000 


000 


768 


108 


903 


2 


916 


617 


106 


250 


236 


2 


916 


617 


677 


766 


815 


1 


636 


801 


016 


875 


796 


25 


000 


000 


010 


704 


806 


-0 


000 


002 


561 


041 


945 


3 


052 


460 


181 


107 


163 


3 


052 


458 


493 


794 


680 


1 


636 


802 


260 


485 


613 


27 


000 


000 


121 


766 


810 





000 


304 


331 


105 


631 


3 


182 


573 


077 


654 


419 


3 


182 


500 


734 


139 


942 


1 


636 


838 


163 


687 


941 


29 


000 


001 


160 


891 


170 


-0 


000 


115 


325 


875 


665 


3 


307 


270 


472 


895 


Oil 


3 


307 


429 


480 


650 


727 


1 


636 


722 


556 


978 


399 


31 


000 


002 


459 


204 


911 





000 


453 


485 


791 


883 


3 


427 


756 


436 


953 


774 


3 


427 


804 


635 


093 


639 


1 


636 


775 


103 


185 


396 


32 


999 


974 


115 


691 


302 


-0 


000 


709 


199 


538 


930 


3 


548 


166 


414 


577 


438 


3 


544 


090 


785 


080 


888 


1 


638 


653 


874 


804 


282 


35 


000 


074 


767 


878 


694 


-0 


001 


304 


748 


848 


729 


3 


658 


326 


733 


913 


608 


3 


656 


678 


487 


565 


035 


1 


637 


404 


037 


276 


759 


37 


000 


475 


686 


981 


808 


-0 


003 


970 


214 


278 


542 


3 


782 


104 


880 


243 


939 


3 


765 


899 


804 


871 


662 


1 


643 


074 


413 


780 


297 


39 


002 


630 


249 


149 


078 


-0 


058 


636 


365 


021 


995 


3 


918 


894 


728 


959 


541 


3 


872 


039 


883 


745 


950 


1 


653 


311 


492 


871 


477 



TABLE II. Numerically calculated eigenenergies E and expectation values (dH/da), (Rcom), and (r re i) of singlet states with 
even parity and vanishing angular momentum for uj = 1, k = 0, and 7 = 1. For this special choice of the parameters uj, k, and 
7 all odd integer values represent exact eigenenergies and the radial expectation value equals (r re i) = 2(2 + v2~7r)/(3 + v27r) = 
1.636 801 341 900 272 (rounded to the last digit). 



Lowest states: 



spin 


parity 


(E) 


{dH/da) 


(-Rcom) 


(rrei) 


(cos ^12) 


singlet 
singlet 
triplet 
triplet 


even 
odd 

even 
odd 


3.445 300 210 741 99 
8.374 386 408 855 55 
5.745 411 335 545 50 
5.526 039 089 152 97 


< IO -13 

< 10~ 12 

< 10~ 13 

< 10~ 13 


0.547 213 052 020 00 
0.953 968 643 580 06 
0.778 149 902 559 62 
0.802 685 965 113 66 


1.410 760 235 037 41 
1.976 361 508 431 29 
1.673 248 129 799 30 
1.722 064 561 651 52 


-0.217 712 873 312 13 
-0.033 320 618 273 03 
-0.066 974 875 946 45 
-0.064 598 448 422 64 


50th excited states: 


spin 


parity 


(E) 


(dH/da) 


(-Rcom) 


(r r el) 


(cos ^12) 


singlet 
singlet 
triplet 
triplet 


even 
odd 

even 
odd 


23.889 095 236 563 
31.288 073 487 
27.856 833 862 282 
27.269 244 649 677 


< io- 12 

< io- 11 

< 10" 12 

< IO -12 


1.357 034 558 
1.494 635 
1.246 927 703 940 
1.424 229 868 91 


2.799 114 48 
3.035 1 
2.427 175 802 70 
2.915 846 906 8 


-0.034 186 244 183 
-0.015 942 
0.026 987 470 565 9 
-0.024 623 600 048 



TABLE III. Zero angular momentum eigenenergy E and expectation values (dH/da), (Rcom), (r r ei)> and (costp^) of the 
lowest state (upper part of the table) and of the 50th excited state (lower part of the table) in each symmetry class. The 
expectation values were calculated from the operator representation by the same method as the Hamiltonian was constructed. 
We round by the last digit, which coincides for at least two different values of a and two different values of nb ase . 



harmonic frequency uj, non- integer values of 7 and an ex- 
ternal magnetic field perpendicular to the plane applied 
to the quantum dot 31 . 

For example, for singlet spin symmetry, even parity 
and vanishing angular momentum the analytical expres- 
sion for one of the energies of the relative motion is 
E rc \ = 2 and the associated radial expectation value is 
(r re i) = 2(2 + v / 27r)/(3 + V^n)- In combination with 
the solutions of the COM motion the total energy reads 
E = E Te \ + Eqom = (2 n + 1), with n € N. We compare 



our results for an optimal choice of a = 0.2 and a ba- 
sis size of n = 6370. For the best converged, low-lying 
values we obtain results in accordance with the analyti- 
cal results up to numerical accuracy of 15 digits. For the 
worst converged values fulfilling our criteria |^S Q | < 0.1 
we still obtain four valuable digits for the eigenenergies 
and two for the radial expectation values (see Table II). 
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2. Full potential case 

We consider the full Hamiltonian (1) with parame- 
ters u> = 1, k — 0.1 and 7 = 1. In this case there 
are no analytic solutions. The convergence of zero an- 
gular momentum eigenvalues is analyzed through their 
dependence on the parameters a £ {0.1,0.15,0.2,0.25} 
and n baso £ {80,90,100,110,120,130}. The basis size 
ranges from 2560 to 13379. For each of these parame- 
ters the eigenvalues and the associated expectation values 
(dH/da), (^com), ('"Vei) an d (cos V12) have been calcu- 
lated. We use a coincidence criterion to establish the con- 
vergence of these quantities: these are converged if they 
are obtained for at least two different values of a and two 
different values of ribaso- The number of coincident sig- 
nificant digits in this process provides the accuracy of the 
result. The lowest eigenvalue and the 50th excited state 
of the four symmetry classes and their respective expec- 
tation values are summarized in Table III. In this numer- 
ical experiment, eigenvalues satisfying (\dH/da\) < 0.1 
typically exhibit at least three converged digits. 



III. ROLE OF THE ANHARMONICITY FOR 
THE COMPLEXITY OF THE SYSTEM 

A. Classical effects 

It is a remarkable property of the system (1), that only 
the case with k ^ and 7 7^ leads to a chaotic classi- 
cal dynamics, which stems from the fact that the system 
separates in two different ways. It separates into a sys- 
tem of two independent particles as long as 7 = and 
into COM and relative motion as long as k = 0. In both 
cases the angular momentum and the energy in the sub- 
systems are preserved. The system with originally four 
degrees of freedom therefore has four constants of motion 
and the underlying dynamics is integrable. Consequently 
we will focus on the full potential case in our classical 
analysis and show some signatures of the chaotic dynam- 
ics. The total energy and the total angular momentum 
remain constants of motion and the phase space is effec- 
tively a six-dimensional space. Except for special cases, 
where the motion is further confined by the choice of the 
initial conditions, the dynamics cannot reasonably be vi- 
sualized by Poincare's surfaces of section. Alternatively, 
insight into the complexity of the system can be gained 
by an appropriate analysis of the dominant frequencies 
of the trajectories 32 . Such frequencies can be identified 
with the frequencies related to the largest weight in the 
Fourier transform of properly selected dynamical quan- 
tities. Associated to each degree of freedom there is 
a dominant frequency, which coincides with the funda- 
mental frequency for regular dynamics. Therefore, com- 
pared with the method relying on the Poincare surfaces 
of section, the frequency analysis is more appropriate for 
highly dimensional systems. Its applications include in- 
vestigations on the Stark-Quadratic-Zeeman problem 33 , 



the two-dimensional standard map 34 , different multidi- 
mensional systems 32 , and the study of the stability of 
the solar system 35 . 

The classical Hamilton function in COM and relative 
coordinates is our basis for the numerical calculation of 
trajectories. The equations of motion are integrated with 
the widely used leapfrog method and the convergence of 
our results is tested with a fourth order Runge-Kutta- 
Nystrom algorithm 36 . We consider trajectories in the 
time interval [t , t + T] with t = 100 and T = 400. The 
energy and the angular momentum is conserved with a 
relative error of the order 10 -5 or better. For each de- 
gree of freedom we consider the combination Xj + ipj of 
the coordinate and momentum, respectively. This quan- 
tity is multiplied by the widely used Hanning Filter 32 to 
avoid effects at the edges of the time interval. A discrete 
Fourier analysis on the above given time interval provides 
the frequencies with an accuracy A/ = 1/T limited by 
the finite time interval. The exact frequency value at the 
peak is found by an iterative scheme in the vicinity of the 
highest values of the discrete Fourier coefficient (golden 
section search). In each iteration the Fourier transform 
is obtained by straightforward quadrature of the Fourier 
integral 32 . 

For our analysis we choose initial conditions such that 
one electron starts at rest from a point on the x-axis, 
while the other electron starts from equally distributed 
positions on the circle with radius one and momentum 
pointing outwards. Then the total angular momentum 
naturally vanishes and we investigate different energy 
regimes. 

The regularity of the system of interacting particles in 
the low energy regime is characterized by the smooth be- 
havior of the dominant frequencies with respect to the 
initial conditions (upper panel of Fig. 1). In this case 
the center of mass is subject to approximately harmonic 
oscillations and the associated frequencies (squares and 
crosses) are constant and mostly degenerate. The fre- 
quencies of the relative motion (circles and pluses) change 
only little or can even be constant in some intervals. The 
dynamics is thus confined to regular nearly harmonic is- 
lands. 

With increasing importance of the anharmonic poten- 
tial for higher energies the dynamics is mixed regular- 
chaotic. This is intuitively clear from the symmetry prop- 
erties of the system and can be verified in terms of the 
analysis of the fundamental frequencies and their sen- 
sitivity to initial conditions. A typical scenario for the 
frequencies in this case is shown in the lower panel of Fig. 
1 . An increase in energy leads to a discontinuous behav- 
ior of the frequencies and results in a complete lifting of 
the degeneracy of the two frequencies of the COM mo- 
tion, which is a consequence of the coupling of the COM 
and the relative motion. 
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FIG. 1. (Color online) Fundamental frequencies of the planar 
classical dot with u = 1, 7 = 1 and k = 0.1 and vanish- 
ing angular momentum. Depicted are the frequencies asso- 
ciated to the relative motion in x- (circles) and y- direction 
(pluses), and to the COM motion in x- (squares) and y- di- 
rection (crosses), respectively. A clearly regular motion is 
recognized in the low energy regime (upper panel: E = 3), 
while for the higher energy a mixed regular-chaotic dynam- 
ics (lower panel: E = 15) is observed. The initial conditions 
are chosen, such that one particle starts static from xi = 0.5 
and yi = 0, while the second starts from equally distributed 
positions on a circle of radius r? = 1 at the angle ip2 with 
momentum pointing outwards with varying absolute value to 
compensate the difference in Coulomb energy. 



B. Quantum effects 

Neither the combination of Coulomb interaction and 
harmonic potential, nor the combination of harmonic 
and anharmonic potential, led to chaotic classical dy- 
namics. In the former case the relative and COM mo- 
tions are separable, while in the latter case the two 
particles are independent. Only a coupling of these 
motions produces a significant impact on the complex- 
ity of the anharmonic two-electron quantum dot prob- 
lem. In this section we investigate the effect of the 
anharmonicity on the complexity of the quantum sys- 
tem. For that purpose we rely on tools that are all 
related to the universal predictions of random matrix 
theory (RMT). These include the next-neighbor-spacing 
distribution P{s) and the distribution P ac {c) of the en- 
ergy gaps c of the avoided crossings which appear by 
slow variation of the parameter k. While the univer- 
sality of RMT predictions has been confirmed by sev- 
eral experimental 3 ' , semiclassical and numerical results 
for systems with underlying chaotic classical dynamics 38 , 
its implications for systems with mixed regular-chaotic 



phase space - such as our system - are still subject of re- 
search and not entirely understood. In the latter case, the 
next-neighbor-spacing distribution P(s) can, with some 
exceptions 39 , be better described by the phenomenolog- 
ical Brody-distribution 40 than by the physically moti- 
vated Berry- Robnik- distribution 41 . A remedy is found 
by extending the idea underlying the Berry-Robnik dis- 
tribution, which is to split the classical phase space into 
distinct regular and chaotic regimes. The effects of dy- 
namical tunneling 42 and chaos assisted tunneling 43 con- 
nect these two classically distinct regimes in the quantum 
regime and already improve the results 44 ' 45 . Recent addi- 
tional achievements considering the effects of flooding 46 
appear to complete the discussion, at least for the next- 
neighbor spacing distribution. Nevertheless a very thor- 
ough analysis of the classical phase space is necessary to 
obtain an ab initio description of this distribution. To 
our knowledge, this has only been performed for one- 
dimensional systems. The distribution P a c(c) is, in the 
mixed case, the sum of a 5-peak, representing real cross- 
ings in the regular regime, and a normal distribution for 
the avoided crossings in the chaotic regime 47 . 
The subtlety of the symmetry and their corresponding 
ensembles, which is fundamental for random matrix the- 
ory, will not be further discussed here. The system under 
consideration belongs to the Gaussian orthogonal ensem- 
ble, as it is symmetric under time-reversal and under ro- 
tations, furthermore the Hamiltonian is real symmetric. 
Our analysis will focus on the dependence on the anhar- 
monicity of P{s), Pac(c) and a measure of the localiza- 
tion of eigenstatcs. The results presented in the following 
were derived from data collected in different runs with the 
numerical method described above. We limit our calcu- 
lations to zero angular momentum for all four symmetry 
classes (singlet/triplet spin symmetry and even/odd par- 
ity). The basis size was determined by the parameter 
"base = 130, which leads to basis sizes of approximately 
12 000 basis vectors, where exact numbers depend on the 
symmetry class. Each calculation supplied more than 
1000 well converged eigenvalues. The quartic potential 
strength was in general varied from k = 0.0 to K = 0.1 
in steps of An = 10 -2 . For the detection of the avoided 
crossings the step size was decreased to Ak = 10~ 5 to 
enhance the resolution of very narrow avoided crossings. 
The universal predictions of random matrix theory can 
only be confirmed if the characteristic spectral proper- 
ties of a system are brought to some general footing. This 
procedure is called "unfolding" and different methods are 
in use. The common two steps for all methods are the 
following: (i) The level density is smoothed by a simple 
fit or semiclassical analysis and (ii) a new set of ener- 
gies is derived from the smoothed level density, in a way, 
that the mean level density is normalized to unity. We 
fit the cumulated level density by a cubic polynomial (i) 
and take the value of the cumulated smooth level den- 
sity evaluated at the former eigenenergy to be the new 
energy (ii). Alternative methods 38 produce qualitatively 
identical results. 
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1. Nearest-neighbor-distribution 

The universality of chaotic properties in quantum me- 
chanical systems has been shown for many examples in- 
volving the statistics of the separations of neighboring 
energy levels. Once the limits for regular (Poisson dis- 
tribution) and purely chaotic (Wigner distribution) were 
observed in different experimental and numerical studies, 
the interest in situations with a mixed classical phase 
space arose and tools to describe the smooth transi- 
tion between these regimes were developed. If there are 
any symmetries left in the problem, the next-neighbor- 
distancc distribution will most likely be a Poisson distri- 
bution. Otherwise, the level statistics of an underlying 
mixed phase space will exhibit a level repulsion which, 
however, is not as pronounced as in the pure chaotic 
case. Alternatives for modeling these statistics are the 
Brody-distribution 40 characterized by the parameter /3, 

PBrody(s) =(0+ l)asP exp(-a^ +1 ) , 



/3+1 



with a = r ! the Berry-Robnik distribution 41 

characterized by the parameter pbr, 



(1 - Pbr) 2 erfc(— ^pbrs)+ 



Pbr(s) = 

+ ( 2/0br(1 - Pbr 

x e^TPlR 



2 

(I-Pbr)s . 



Pbr- 



and the extension of the latter proposed by Podolskiy 
and Narimanov 45 ' 48 characterized by the parameters p 
and V° c , 



Ppn(s) 



(W) 'F[ 



2p{\-p)F 



V RC 



+ g P 3 I x 



xe *' 



-(l-p)s 



(16) 



with 



F(x) 



1 



The Brody parameter (3 describes the transition from 
regular (f3 = 0) to chaotic (j3 — 1) behavior. It must 
be noted, that the Brody parameter j3 lacks a quan- 
titative physical meaning, but describes the transition 
merely qualitatively. The Berry-Robnik distribution is 
characterized by the parameter pbr, which is the ra- 
tio of the chaotic to the total phase space volume and 
is a purely classical property of the system. The gen- 
eralization of the Berry-Robnik distribution by Podol- 
sky and Narimanov include perturbativc quantum cor- 
rections through the parameter V^ c , which describes 




k [ a . u .] 



FIG. 2. (Color online) Left: The numerical level spacing dis- 
tribution (gray bars) for singlet states of odd parity is shown 
together with the fitted Brody distribution (solid line), the fit- 
ted Berry-Robnik distribution (dot-dashed line) and the fitted 
Podolskiy-Narimanov distribution (dashed line) for the case 
with k = 0.02. Right: Brody parameter f3 (circles), Berry- 
Robnik parameter pbr (squares) and Podolskiy-Narimanov 
parameter p (crosses) singlet states of odd parity for several 
values of k 6 [0, 0.1]. 



the tunneling between regular and chaotic regions and 
chaos assisted tunneling between regular regions via the 
chaotic sea. 

We have considered all states from the principal quan- 
tum numbers n = 1 through to n = 40, in the unper- 
turbed system, for all symmetries, leading to a specific 
number of eigenenergies according to Table I. In Fig. 2 
we show the results of the statistical analysis for the sin- 
glet spin symmetry odd parity case performed for 715 
eigenvalues (a = 0.2, ribaso = 130, n to t — 11 168). In the 
left panel our numerical data, for the specific case with 
k = 0.02, is shown as gray bars, while the different fits 
according to the previously described distributions are 
given by the solid, dashed and dot-dashed lines. For the 
Berry-Robnik distribution we fitted the parameter pbr- 
The fit for the numerical data using the Brody distribu- 
tion is in general better than the fit by the Berry-Robnik 
distribution, but has the disadvantage of an unphysical fit 
parameter. The Podolskiy-Narimanov distribution im- 
proves the quality of the fits, naturally, as a second fit 
parameter V^ c is introduced. For a quantum dot with in- 
teracting electrons (7 = 1) without quartic term (n = 0) 
we analyzed the spectrum of the radial equation of the 
relative motion. Being still too close to the pathological 
harmonic oscillator case no proper fit could be performed. 

For the full potential the Brody-parameter acquires an 
intermediate value of /3 « 0.45 (/3 sa 0.2) for the singlet 
even parity case (other cases). Variations of the strength 
k of the quartic potential leads for k 6 [0.01,0.1] only 
to small variations of /3 (see right panel of Fig. 2). The 
parameter pbr exhibits a similar behavior however its 
value is shifted with respect to f3 by an approximately 
constant value 0.3. This stems from the similarity of 
the two distributions with this special shift in the pa- 
rameters. The parameter p of the Podolskiy-Narimanov 
distribution tends to be close to one of the two previous 
parameters. The quantum mechanical coupling parame- 
ter Vft C lies between 0.1 and 0.45. 
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Simple algebraic considerations show that the strength 
of the quartic potential only shifts the energy regimes for 
which the anharmonic effects become important. There- 
fore, as long as the analysis includes a regime wide 
enough in energy the effect of the anharmonicity can be 
measured for any non vanishing value of k and should not 
depend strongly on it. Performing the statistical analy- 
sis for a lower number of eigenenergies, cutting off at a 
specific principal quantum number n < 40 leads to a de- 
crease of the chaoticity parameters for k < 0.05, while 
the results for k > 0.05 remain unaltered. 



2. Avoided crossings 

A manifestation of the non-separability of a system and 
its complexity is the occurrence of avoided crossings in 
the spectrum depending on a slowly varying parameter. 
This is the case in our system when k is the adiabatic 
parameter (see upper panel in Fig. 3). Avoided cross- 
ings are naturally related to a drastic interchange of the 
symmetry properties of the eigenstates involved. This is 
typically characterized by abrupt changes in the behavior 
of some expectation values. For instance, as the eigen- 
states \n c , m c , n r , m r ) ep = |l,4,0,-4) + and |0,2,3,-2) + 
of the harmonic (k = 0) quantum dot evolve as n varies, 
there is a sudden jump in the expectation value of the 
COM radial distance (R) close to the avoided crossing 
around k = 0.00716 as shown in the lower panel in Fig. 
3. Less pronounced is the change of (R) for the states 
|0,0, 5,0) + and |3, 2,0,-2) + close to the avoided cross- 
ing around k = 0.00325, though the respective wave 
functions 49 completely interchange their properties as 
seen in the side insets of Fig. 3. 

For the construction of the distribution P a c(c) we re- 
quire an efficient determination of avoided crossings in 
a large amount of spectral data. Though jumps in the 
expectation values (R) can be used to identify avoided 
crossings, the systematic detection of these jumps is not 
a trivial task. Alternatively, the quantum fidelity sus- 
ceptibility of the eigenstates provides an efficient method 
for such purpose. The quantum fidelity susceptibility 
X of an eigenstate ip n is equivalent to the curvature of 
tp n depending on the varying parameter 00 . It can be 
calculated via the static quantum fidelity Fs K (K,n) = 



\{ipn,K I i>n,K+6K.)\ , 



X 



lim 

<5k^0 



i-f Sk 

(5k)* 



lim 

<5k->0 



log(Jfa) 
(6k)* ' 



and is largely independent of the perturbation i5k 51 . 

The typical behavior of the quantum fidelity suscep- 
tibility close to an avoided crossing is illustrated in the 
middle panel of Fig. 3. It is characterized by three prop- 
erties which are easy to implement for practical purposes: 
The susceptibility has a peak near an avoided crossing, 
this peak is nearly identical for the two non-crossing 
states and the mean values of the susceptibilities before 
and after the peak interchange for these two states. 
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FIG. 3. (Color online) Upper panel: Example of two avoided 
crossings in the regime of weak quartic potential. The plot 
shows eigenenergies for four singlet states with even parity 
and vanishing total angular momentum in the potential with 
uj = 1.0, 7 = 1.0 and varying k. The eigenstates are la- 
beled by the exact quantum numbers n c , m c , n r , and m r of 
the harmonic problem (k = 0). Middle panel: the quantum 
fidelity susceptibility for the same states. The susceptibili- 
ties of the two non-crossing states show peaks, which coincide 
at the maximum and are interchanged after the peak. Lower 
panel: The expectation values of the center of mass radial dis- 
tance (R) in the vicinity of the avoided crossings of the upper 
panel. Side panel: Contour plots of the states |0, 0, 5, 0} + and 
|3, 2, 0, -2} + for k = 0.002 before the avoided crossing ((a) and 
(c), respectively), and for k — 0.005 after the avoided crossing 
((b) and (d), respectively). 



After checking that the susceptibility is indepen- 
dent of the perturbation for several values of 6k <G 
{10~ 6 , 10~ 7 , 10~ 8 , 10~ 9 } data collection was performed 
using 5k = 10 -7 . With this method it was possible to 
determine nearly 9000 avoided crossings within a range 
of 0.0 < k < 0.01 and energies up to 70. The levels 
start to mix depending on the quartic potential strength 
k, though for small values of k the mixing can only be 
observed for high lying energy states. 
For each avoided crossing the energy gap between the 
two close lying states, known as the width of the avoided 
crossing, is calculated and a statistical analysis of these 
values is performed. The expected distribution is the 
weighted sum of the distribution of the widths of the 
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FIG. 4. Fitting parameter A for the distribution of the widths 
of the avoided crossings for approximately the first 1000 eigen- 
states of the system. The right plot shows an example for the 
numerically determined distribution (grey bars) in the region 
0.002 < k < 0.003 and the fit Eq. (17) without the <5-peak 
close to zero (black line). 



avoided crossings for the chaotic case with a (5-peak, 
which represents the non-avoided crossings of the regular 



cas(> 



17. 



P(c) = (1 - X)S(c) + 



2A 2 



■ exp 



-A 2 c 2 



(17) 



where (c) represents the mean value of the widths of the 
avoided crossings. We fit with the cumulated distribution 
following Rcf. 47. 

Even for small values of the quartic potential strength 
a high value for the fitting parameter A indicates a high 
degree of chaoticity in the system, see Fig. 4. Again the 
very onset of the anharmonic term seems to carry the 
main effect, while the value of the interaction strength 
does not significantly change the distribution. The same 
argument as in the previous section holds here, as we an- 
alyze a wide energy regime capable to represent the fea- 
tures of anharmonicity even for small values of k. The 
number of avoided crossings is already high for small val- 
ues of K, because the levels, originally degenerate in the 
purely harmonic case, strongly mix as soon as k ^ 0. 



3. Eigenvector localization 

We calculate an eigenstate specific information entropy 
depending on the parameters k and 7, which depends on 
the basis representation of the eigenvectors. Following 
the arguments by Zelevinsky 52 we obtain physically rel- 
evant results by choosing a reference basis for the rep- 
resentation of the eigenvectors that is physically related 
to the system under consideration. In our case this is 
the purely harmonic oscillator described in section II B. 
A general eigenvector of this system is given by 



\lpnh 



G {0,..,N},ke{l,..,K n }, 



where N goes to infinity and the values of K n are given in 
Table I. We calculated numerically the harmonic case in 
order to achieve a complete (K n ) basis representation for 



all results presented here. Instead of determining the en- 
tropy directly in the harmonic oscillator basis, we project 
on the energy subspaces. For a general vector \ip) , it holds 

00 K n 00 

\<p) = ^2^2c nk \if; nk ) = 2_,e n \E n ), 

n=0 fe=l n=0 



with 



\E n ) = — y^Cnfcl^nfc) 



k=l 



and the coefficients 



\ 



K, 

E 

fc=i 



\C n k\ 



chosen such that the representation basis \E n ) is or- 
thonormal. 

The information entropy defined by 



JV 



^ = 5> n | 2 io g (|e„| 2 ), 



is a measure of the localization of the vector \ip) in the 
harmonic oscillator basis. Large values of S v imply a 
large spread of the state \ip) in this basis. 

Fig. 5 shows the entropies of cigenstates with singlet 
symmetry and even parity of a two-electron quantum dot 
with a harmonic confinement (ui = 1) and four different 
situations. As our reference basis is the harmonic oscil- 
lator basis the purely harmonic case gave only vanishing 
entropies and is not considered here. For the harmonic 
case with electron-electron repulsion, in (a), a very regu- 
lar behavior of the entropies can be observed. Indeed in 
this case the entropy is only a measure for the relative 
motion and the horizontally ordered entropies belong to 
states which differ only in the COM quantum numbers. 
The eigenstates contained in Table II are highlighted with 
circles and systematically belong to the states within the 
region of higher entropy. These states are configurations 
with no angular momenta in the subsystems of the COM 
and of the relative motion. The second class of states, 
showing lower values for the entropy have non- vanishing 
angular momenta m c and m r = —m c in the subsystems. 
We have marked a series of these states with diamonds. 

In Fig. 5 (b) we show the entropies with weak anhar- 
monic confinement (k = 0.01) and Coulomb repulsion. 
The behavior of the lowest eigenstates doesn't change 
significantly, up to an energy of£w 10, while above this 
limit the entropies increase and lose most of their regular 
structure. The former effect can clearly be related to the 
narrowing of the potential and the higher number of har- 
monic states necessary to represent the eigenstates. The 
latter is a signature of the rising complexity of the system, 
induced by the coupling of the two previously separate 
motions and the increasing number of avoided crossings. 
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FIG. 5. (Color online) Information entropies of states with 
singlet symmetry and even parity for four combinations of 
the parameters re and 7: (a) re = 0.1, 7 = 0, (b) re = 0.1, 
7 = 1, (c) re = 0.01, 7=1, and (d) re = 0, 7 = 1. In 
all cases there is a harmonic confinement with ui — 1. The 
entropy has been calculated for approximately 1200 eigen- 
states in (a), 700 eigenstates in (b) and 250 eigenstates in 
(c) and (d). In (a) we marked the states with the quantum 
numbers \n c ,m c ,n r ,m r } ep — |fc,0,0,0} + with circles (also in 
(b)), \k, 0, 5, 0) + with squares and |fc, 2, 1, -2) + with diamonds, 
where k is an integer number. The other symmetry classes 
exhibit similar behavior. 



The occurrence of avoided crossings can in particular be 
seen from the jumps in the entropies of the states marked 
with circles, which we obtained by adiabatically following 
the levels from (a) to (b). 

These effects arc further enhanced for the stronger an- 
harmonic (re = 0.1) case in (c), where only the ground 
state entropy is close to the harmonic case. This is in con- 



trast to the case (d) with anharmonic confinement but 
without electron-electron interaction. The values of the 
entropies are as high as in case (c), but they still show 
a specific structure. This can be understood by consid- 
ering that, in this case, the system is again separable in 
two independent particles. 

Our analysis of the information entropy shows again, 
that the complexity stems from the interplay of the inter- 
particle interaction and the anharmonic confinement. As 
long as one of those is omitted a high degree of regularity 
is still present, reflected in the behavior of the entropies. 
Qualitatively similar observations are obtained when we 
choose for instance the basis (10) as a reference. 



IV. CONCLUSION 

We have presented a detailed description of a numer- 
ically exact treatment of planar two-electron quantum 
dots, which extends the common harmonic model by in- 
troducing a quartic potential term. Our approach pro- 
vides an accurate characterization of the spectrum of this 
system for a wide range of different parameters and sym- 
metry classes. This has been exploited for studying an- 
harmonic effects in the complexity of this system. Our 
analysis showed that generally the interplay between the 
Coulomb interaction and the anharmonic term is respon- 
sible for a significant reduction of the regions of regular 
classical motion of the planar quantum dot, which oth- 
erwise is integrable. On the quantum level, signatures of 
this mixed regular-chaotic underlying dynamics are ob- 
served in the level spacing distributions, the appearance 
of avoided crossings and the distribution of their separa- 
tions, and in the eigenvector localization entropies. The 
complexity arises as soon as the anharmonicity is intro- 
duced in the interacting harmonic quantum dot. The 
complex quantum regime at rather high energies in the 
case of small anharmonic perturbations is shifted to lower 
energies as the anharmonicity increases. Apart from this, 
a common feature observed is the rather low dependence 
of the studied complexity measures on the strength re of 
the quartic confinement. 

Our approach can be readily used for the study of fur- 
ther phenomena in two-electron quantum dots including 
the effects of the anharmonicity in Hund's rules and the 
consequences of chaos in the entanglement of two elec- 
trons. Our approach is also suitable for the study of 
decoherence processes in quantum dots. For instance, 
the anomalous behavior of quantum fidelity decay found 
in many-body systems 53,54 can be tested for a complex 
anharmonic system without assuming mean- field approx- 
imations. Furthermore, the algebraic representation of 
observables, in particular the dipole operator, can be eas- 
ily implemented. This offers the opportunity to investi- 
gate phenomena connected to the interactions of quan- 
tum dots with laser pulses. 
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Appendix A: Modified atomic units 

We apply modified atomic units (7 = 1), where we set 
h = to* = (47reoe) = q e = fee = 1- The natural scales 
are the modified Bohr-radius 

4:ire eh 2 
o-q = 5—, 

the modified Hartrcc energy 



the natural frequency 

_ Eh_ _ rn*q* 
VQ ~ h ~ h 3 (47re e) 2 ' 

and the natural timescale 

h _h 3 (47re e) 2 

*0 — ~5 ~ 3 • 

Eh m*ql 

Solid state quantum dots are most commonly realized 
on substrates of Indium Arsenide (InAs) and Gallium 
Arsenide (GaAs). For both cases we supply the values of 
the effective electron mass m„, the dielectric constant e 
and the natural scales: 

• InAs: to* ps 0.023 m e , e ps 15.15, ao ps 35 ran. 
E h ps 2.7 mcV, u ps 4.14 THz, £ « 0.24 ps, 

• GaAs: m* ps 0.063 m e , e ps 12.9, do ps 11 nm, 
£ h ps lOmeV, i/ ps 15.65 THz, i ps 0.066 ps. 

A typical value for the confining potential is 
Eq = hajQ = 3 meV , Ref. 22, which implies that uj 
is of order unity. The ratio of the harmonic confining 
potential to the natural energy scale of one Hartree 
(i?h) is expressed as u> = E^/E^. We set oj = 1 in our 
calculations. 



Appendix B: Coordinate basis representation 

In the polar coordinates associated to the parabolic coordinates (5) the basis functions can be represented in 
coordinate space with identification of the quantum numbers: 

M = ni — n 2 j L = 713 — 714, N = min (m, 712) , A' = min (713, 714) , 

as follows: 

<p n (r + ,4> +> r-,<t>-) = (m,n 2 ,n 3 ,n4 | r+, 0+,r_, 0_) = A^V* 1 iif (r 2 ) (r 2 ) e -iW-^_) e -i W++ i*-) ) 
where are associated Laguerre polynomials and j\f is a normalization constant: 



N = {-l) 



N+K 



Nl K\ 



it V (|M|+JV)! (|L| + A)!' 



The density plots shown in Fig. 3 are calculated by integrating over the angles in COM and relative coordinates, where 
the appropriate coordinate transformation is r\ = 4i?coM and = 2r Te \. For a general vector \ip) = ^2 n a n \n) we 
plot the following function: 

7T 7T 

^2^2 a n a m J d(j) + J d(j>- i?COM ^rel <Pn(V ^ RcOM, </>+, \/4 r rc i, (/)_)* p m ( V^RcOM, </>+, \/4 r rc i, 
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